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ABSTRACT 

Collapsars are fast-spinning, massive stars, whose core collapse liberates an energy, that can be 
channeled in the form of ultrarelativistic jets. These jets transport the energy from the collapsed core 
to large distances, where it is dissipated in the form of long-duration gamma-ray bursts. In this paper 
we study the dynamics of ultrarelativistic jets produced in collapsars. Also we extrapolate our results 
to infer the angular energy distribution of the produced outflows in the afterglow phase. Our main 
focus is to look for global energetical properties which can be imprinted by the different structure of 
different progenitor stars. Thus, we employ a number of pre-supernova, stellar models (with distinct 
masses and metallicities), and inject in all of them jets with fixed initial conditions. We assume that 
at the injection nozzle, the jet is mildly relativistic (Lorentz factor ^ 5), has a finite half-opening 
angle (5°), and carries a power of lO^^ergs"^. In all cases, well collimated jets propagate through 
the progenitor, blowing a high pressure and high temperature cocoon. These jets arrive intact to the 
stellar surface and break out of it. A large Lorentz factor region F > 100 develops well before the 
jet reaches the surface of the star, in the unshocked part of the beam, located between the injection 
nozzle and the first recoUimation shock. These high values of F are possible because the finite opening 
angle of the jet allows for free expansion towards the radial direction. We find a strong correlation 
between the angular energy distribution of the jet, after its eruption from the progenitor surface, and 
the mass of the progenitors. The angular energy distribution of the jets from light progenitor models 
is steeper than that of the jets injected in more massive progenitor stars. This trend is also imprinted 
in the angular distribution of isotropic equivalent energy. 

Subject headings: hydrodynamics - jet - GRBs - supernovae - shock - relativity 



1. INTRODUCTION 

Recent observations of gamma-ray bursts (GRBs) sug- 
gest that long duration GRBs and type Ib/c super- 
nova (SN) explosions are tightly connected. For ex- 
ample, SN1998bw was observed in the pos itional er- 
ror box of GRB980425 ijGalama et al.lll998[ ). In this 
case the GRB/SN association was based on the spatial 
and temporal coincidence of both events. The most re- 
markable example of long GRB/SN link came in 2003, 
when the spectra of both the GRBG3G329 afterglow and 
of the SN2003dh were measured, since the burst hap- 
pened closeby and it was quite bright. The supernova 
spectrum, which includes many complex lines, gradu- 
ally appeared from the decaying afterglow spectrum af- 
ter a few tens of days from the burst. The spectrum 
of SN2003dh after about a month from the explosion is 
quite similar with that of SN1998bw at the same stage. 
Both SN1998bw and SN2003dh are typ e Ic supernovae 
(jlwamoto et al.lll998HStanek et al.ll2003[ ) whose progeni- 
tor had lost the hydrogen and the helium envelopes dur- 
ing the pre-supernova stage. They are also categorized 
within a special class of supernova explosions, so-called 
hypcrnovae, whose explosion energy is about ten times 
higher, i.e., ~ 10^^ erg, than that of ordinary super- 
nova. Indeed, our common view is that most long-lasting 
GRBs are produced by core-collapse supernovae akin to 
SN1998bw. 



Electronic address: 'AM, e-mail: mizuta@cfs.chiba-u.ac.jp' 

^ Center for Frontier Science, Chiba University Yayoi-cho 1-33, 
Inage-ku, Chiba, 263-8522 Japan 

^ Departamento de Astronomia y Astrofisi'ca, Univcrsidad de Va- 
lencia, 46100-Burjassot (Valencia), Spain 



GRBs are not exclusively linked to hypcrnovae. For 
instance, the long lasting burst (i^ ~ 2 x 10^ s) 
GRB060218 (or XRFG6G218) was associated with the 
type Ic SN2GG6aj, which is not a hypernova. The ex- 
plosion energy falls within th e regular range f or typ e Ic 
events (jCampana et al.ll2006D . IMazzali et all (|2GG6f) ar- 
gued that the progenitor star may not form a black hole 
but a neutron star, since the estimated mass of the pro- 
genitor during the main sequence is 2GM0. On the 
other hand, there are recent examples of cases of long- 
lasting GRBs (GRBG60505 and GRBG6G614) where no 
supernova signature was observed at all even if, con- 
sidering the distance, the observati onal trace of a su- 
pernova should have been detected (iFvnbq et al.l [20061 : 
IDella Valle erall[2006l : iGal-Yam et al.ll2006n . 

Though a variety of long duration GRBs which have a 
strong connection with SNs are observed, we still lack of a 
complete picture of the processes by which a sizable frac- 
tion of the energy involved in a type Ib/c SN explosion 
is tapped in a relatively narrow channel, and produces a 
GRB at a large distance from the original site of genera- 
tion. iRees fc MeszarosI (|1992[ ) proposed t hat the d e ath o f 
massive stars can be an origin of GRBs. IWooslevI (IT993h 
introduced the coUapsar model to account for the pro- 
genitor system of long GRBs. According to this model, 
a non-sphcrical outflow could be formed from the deep 
inside of the progenitor where a black hole or a proto- 
neutron s tar is born as a result of the collapse of the 
uon core (jMacFadven fc Wooslevlll999f ). If the specific 
angular momentum of the iron core is sufficiently large, 
an accretion torus may develop around the central com- 
pact object. The system formed by a central object and 
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an accretion disk has the potential to launch a bipolar 
outflow. 

The mechanisms proposed to extract energy out 
of such central engines arc basically two (e.g., 
lAlov fc Obergaulinger.2007i ): thermal or hydromagnetic. 
Thermal mechanisms rely on depositing a considerable 
amount of thermal energy in the vicinity of the rota- 
tion axis of the system, just above the poles of the cen- 
tral compact object, where a low density funnel has 
develop in the course of the evolution. The accre- 
tion energy of the hot torus is converted into a co- 
pious flux of neutrinos v and anti-neutrinos v. From 



the 



-annihilation a hot 



-plasma results. In its 



turn, e"'"e~-pairs annihilate yielding a fireball of high 
energy photons. The conversion of the thermal en- 
ergy of the fireball into kinetic energy partly determines 
the subsequent evolution of the plasma. It accelerates 
to ult rarelativistic spe eds (reaching Lorentz factors of 
~ 50 lAlov et all [2000h while, at the same time, inter- 
acts with the progenitor system. Alternatively, MHD 
process may tap a fraction of the rotational energy 
of the BH or of the accretion disk to form an out- 



flow (IProiaeranJloMM McKinnevI 
20061: iNagataki et al.ll2007t iKomissarov fc Barkovii2007t 



Takiwaki et al.ll2009f ) 



A number of works have dealt with the hydrodynamic 
properties of the outflows generated in coUapsar progen- 
itors as they propagate through the progenitor system 
(and in some cases beyond the surface of the progeni- 
tor star). The problem is addressed by means of nu- 
merical relativistic hydrodyna mic simulations wi t h dif- 
ferent degrees of complexitv JZhang et al.l 120031 [200^: 
Mizuta et a"ni2006tlMorsonv et al.ll2007l: iTominaga et all 



20071) ■ which assume that a quasi-steady momentum flux 
has been produced at a certain distance from the region 
where the energy is released (independent of which is the 
actual energy extraction mechanism -MHD or thermal-) . 
Therefore, such numerical works assume the existence of 
a nozzle through which the injection of a supersonic jet 
is produced, and put their focus on the modification of 
the morphology and of the dynamics of coUimated out- 
flows as they travel through the progenitor star. Depend- 
ing on the exact inflow condit ions, a variety of di fferent 
outflows result. For instance. iMizuta et al.l (|2006f ) finds 
a whole spectrum of outfiows ranging from coUimated, 
relati vistic jets to p oorly coUimated expanding winds. 
Thus, IMizuta et al.l argue that such a variety of result- 
ing outflows supports the idea that the same coUapsar 
scenario can yield a number of different phenomena (in 
agreement with the pr evious ideas of unifi cation of high 
energy transients, e.g.. iRamirez-Ruiz et al . 2002b), such 
as, GRBs, X-ray rich GRBs, X-ray Flashes, and normal 
supernovae. 

In the prompt GRB phase we observe the emission 
of high energy photons from an ultrarelativistic outfiow, 
which is generated at a distance 10^^~^^ cm, very far from 
the central engine. Due to the large optical thickness of 
the outfiow at scales comparable to that of the progeni- 
tor, we have not observed any electromagnetic emission 
directly from the prc igenitor s c i far, except, perhaps, some 
precur sor activity (jWooslevI 120001 : iRamirez-Ruiz et aD 
l2002aD . Since the direct detection of progenitors of GRBs 
is nowadays impossible (if they happen at cosmological 
distances), the observation of the association between 



GRBs and SNe probably provides the best clue to under- 
stand the progenitors of the GRBs. Though more than 
100 GRBs per year are identified, most of them occur 
far away. Thus, it is technically impossible to identify 
their accompanying supernovae. Another clue regarding 
the nature of the GRB progenitors comes from the envi- 
ronment and the host galaxies in which the burst takes 
place. Both clues indicate that long-duration GRBs 
are associated with t he death of the most massive stars 
(jFruchter et al.ll2006f ). The typical hosts of long GRBs 
are star-forming, low metallicity galaxies (with an star 
formation rate ^ 10'^M (7iy~^: iBerger. Kulkarni. fc FraiH 
l200lHFrafl et al.ll2002f ) but bluer than typi cal starburst 
galaxies, with little dust (|Le Floc'hl l2004f l. and lower 
masses than current ellipticals, i.e., they correspond to 
the typical environm ents of formation of massive stars. 
iFruchter et al.l ()2006D also conclude, that the host galax- 
ies of the GRBs are significantly fainter and more irreg- 
ular than the hosts of core-collapse supernovae. 

A very interesting question is weather it is possible to 
get any information on the nature of a GRB progeni- 
tor from the observation of the afterglow emission. One 
possibility is to look at the angular distribution of in- 
tegrated energy per unit of solid angle, as observed in 
the af terglow phase of the burst. iLazzati fc BegelmanI 
(|2005| ) estimated theoretically such an angular distribu- 
tion assuming, that the kinetic energy of the jet is con- 
verted to thermal energy in the cocoon, till the head of 
the jet reaches the progenitor surface. The cocoon orig- 
inates from jet material which crosses through the ter- 
minal strong shock of the coUimated outfiow and moves 
away from the center of the progenitor surrounding the 
beam of the jet. A fat cocoon develops for light jets, i.e., 
jets whose rest-mass density is much lower than that of 
the ambient gas into which the jet propagates. When the 
jet breaks out the progenitor surface, the thermal energy 
is released to a low-density inter stellar medium (ISM). 
ILazzati fc BegelmanI ()2005[ ) concluded, that the energy 
distribution per solid angle {dE/dO) of the jet displays a 
6^^ dependence with the viewing angle af ter its eruption 
throug h the progenitor surface. Recentlv. iMorsonv et al.l 
(|2007| ) have used hydro dynamic simulations to test the 
theoretical prediction of ILa zzati fc BegelmanI . and found 
that their numerical models do not follow the inferred 
theoretical angular energy distribution. 

In this paper we also try to verify the analytic rela- 
tion for the angular dependence of the energy with the 
polar angle that was proposed bv ILazzati fc BegelmanI 
(120051). We explore a pa rameter space different from that 
of IMorsonv et al.l ()2007D in order to compute the depen- 
dence of the angular energy distribution on the structure 
of the progenitor. The progen itor models are bu il t upo n 
the pre-supernova models of IWooslev fc Heged ([20M ). 
Along the way, we also characterize the hydrodynamic 
properties of relativistic jets propagating through differ- 
ent progenitor stars. 

The paper is organized as follows. We describe our 
physical model, the choice of stellar progenitors, and rel- 
evant numerical details in Sect. [21 In the appendix, we 
provide a study of a selected sample of models in order 
to justify our choice of numerical resolution and the ef- 
fects it has on our conclussions. The dynamics of the 
injected bipolar outflows, and the extrapolated angular 
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energy distribution in the afterglow phase is considered 
in Sec. [31 FinaUy, we discuss our results and write down 
the conclusions of this work in Sect. ID 

2. MODEL 

We will investigate the dependence of the properties of 
relativistic jets, injected in a pre-supernova stellar model 
at a certain distance from the center, using relativistic 
hydrodynamic simulations. In Sect. [^TTl we show the dif- 
ferent stellar progenitors used in this study. We pro- 
vide the technical details of the numerical simulations in 
Sect. 12.21 Finally, in Scct. l2.3l wc specify the physical con- 
ditions used to inject relativistic jets in the pre-supernova 
progenitors described above. 

2.1. Progenitors 

In the last years, some detailed calculations of stellar 
evolution of massive stars have been done including the 
effects of initial angul ar momentum, dynamo, metallic- 
ity, and mass loss rate (lYoon fc Lange jl2005l : lYoon et al.1 
mim IWooslev fc Hegedl2'006i r According to these stud- 
ies, the metallicity of the progenitor strongly affects the 
evolution of the angular momentum distribution at the 
pre-supernova stage, in such a way, that low metallicity 
is preferred to obtain a large angular momentum in the 
core of the progenitor. 

For the purposes of this work, we e mploy some of the 

fi-supernova models computed by IWooslev fc Hegerl 
10(3). We stick to the same naming convention than 
the former authors, and consider several sets of models 
(Tab.[T]). The f irst group correspo nds to the IIE16-series 
of 16 models of lWooslev fc Hegeil . which include progen- 
itor stars for which WMq bare helium cores are evolved, 
that have solar metallicity, and different amounts of ini- 
tial angular momentum, dynamo effects, mass-loss rates, 
etc. The last three models of Tab. [H 160C, 16TB, and 
16TC form the second group of progenitors. They cor- 
respond to stars with the same initial mass as those of 
the first group but with a smaller metallicity {Ziqqc = 
O.IZq, Zi6tb = ZieTC = O.OIZq). The second set of 
low-metallicity models has been chosen among the many 
other possibilities available because their radius, at the 
pre-supernova stage, are the smallest among all other 
low-metallicity progenitors (in all cases, their stellar radii 
are i?, < 10^^ cm). 

Figure [T] shows the radial mass profiles of the mod- 
els HE16C, HE16L, and HE16N. HE16C is representa- 
tive of progenitors whose pre-supernova mass is small 
-Mhei6C — SMq (due to the vigorous mass loss rate in 
the late phases of its evolution; Tab. [1]). Other mem- 
bers of this group of low pre-supernova mass arc IIE16B, 
IIE16J, and HE16K. We will refer to this group as type- 
L. The model IIE16N belongs to the group of more mas- 
sive progenitors (Mhei6C — ISM©), to which we will 
refer as typc-H models. Finally, the model IIE16L (also 
IIE16D) falls in the middle of these two groups (mem- 
bers of this group will be called type-M models). Its total 
mass is about 9.5M0. Figure [2] shows radial mass den- 
sity profiles of the low metallicity models 16TB, 16TC, 
and 160C. Though the mass of the models is similar 
(~ I5.3M0), the density profiles are slightly different. In 
total, 19 models are considered in this study (Tab. [T|). 
We neglect for the progenitors any deviation from spher- 
ical symmetry arising from the rotation of the models. 
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Fig. 1. — Radial mass profiles of models HE16C, HE16L and 
HE16N. 
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Fig. 2.— Radial mass profiles of models 16TB, 16TC and 160C. 



Therefore, for the study we present here, each progeni- 
tor differs from the other mainly in its total mass, radius, 
and mass density profile at the pre-supernova stage. 

Though a nonspherical structure is expected around 
the black hole due to the rapid rotation of the progenitor, 
it is reasonable to assume spherical symmetry for the en- 
velopes of the progenitor for radial distances r > 10^ cm, 
which is where we put the innermost radial boundary in 
the numerical simulatio ns of this study. Thus, we only 
take from the models of lWooslev fc Hegerl (|2006f ) the ra- 
dial density and the radial velocity profiles which result 
by the end of the pre-supernova evolution. We assume 
that the pressure of the progenitor is very low, or, equiv- 
alently, that the initial specific internal energy (e) is set 
to be very low (e/c^ = 10^®, where c is the speed of 
fight). 

Both the gravitational force produced by the central 
compact object, and the progenitor self gravity are ig- 
nored, since the timescale for the outfiows to cross the 
progenitor and to break out from the surface of the pro- 
genitor, 3s, is much shorter than free fall timescale for 
the stellar envelopes. 

We extend the radial mass density profile to the outside 
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of the progenitor up to the outer computational bound- 
ary located at r^ax = 10^^ cm. For the models which 
barely loose mass during the latest stages of their evo- 
lution, the rest-mass density is assumed to be uniform 
(piSM = 10~^ g cm^"^) and much smaller than that at the 
progenitor surface. If the progenitor star has a non-null 
mass loss rate (i.e., if t he parameter a > 0.01, according 
to the nomenclature of lWooslev fc Heg"eill2OO60 . we take 
a dependence in radial mass-density profiles from the 
surface of the progenitor (see Fig. [T]). 

2.2. Computational Domain and Basic Equations 

We map the spherically symmetric progenitor models 
of Sect. [221 into a two-dimensional grid in spherical coor- 
dinates (r X 6). We assume that our models arc axial and 
equatorially symmetric and, therefore, specify reflection 
boundary conditions at the polar axis {9 = 0°) and at the 
equator (6* = 90°). The radial grid consists of Nr = 1000 
points, uniformly spaced in log r, which extends from 
Tmin = 10® cm to r^ax = 10^^ Cm. The smallest radial 
grid spacing, besides r,„in, is Armin = 10^ cm, while the 
largest one, besides rmax, is Ar,„ax = 6.5 x 10* cm. The 
resolution we have choosen here represents a trade-off 
between accuracy and feasibility of the numerical simu- 
lations, as we discuss in the Appendix. Free outflow (i.e., 
zero gradient) boundary conditions are set at r = rmin 
and Tmax- The polar grid has Ne = 180 grid points uni- 
formly spaced in the range 0° < 6* < 90°, {A9 = 0.5°). 
We use the same 2D special rela tivistic hydrodynamic 
code of iMizuta et~all (|2004 12006[ ) to perform our simu- 
lations. The code provides 3rd order accuracy in both 
space and time, by applying a PPM intra-cell interpo- 
lation and a TVD-Runge Kutta time integration. For 
the sake of simplicity, we employ an ideal gas equation 
of state {p ~ (7 ^ l)pf) with uniform adiabatic index 
7 = 4/3, where p, p and e are the pressure and the rest- 
mass density, respectively. 

2.3. Jet Injection Conditions 

We assume that a jet has been generated by the central 
engine, and that at a certain distance, quasi-steady injec- 
tion conditions are settled through a well defined circular 
nozzle. Thus, we inject plasma, in the radial direction, 
through the innermost radial boundary at r = rmin in 
a cone of half-opening angle 6j = 5°. The jet injec- 
tion proceeds for a period tinj = 4 s. We parametrize 
the outflowing plasma by assuming that it is hot (we set 
Cj/c^ = 30) and moderately relativistic (the Lorentz fac- 
tor being Fj^o = 5). We adopt the convention that the pa- 
rameters of the outflow at the injection point are named 
with a subscript 'j'. Because of the conversion of thermal- 
to-kinctic energy, the injected flows have the potential 
to accelerate to bu lk Lorentz factors larger than 100 
(|Mizuta et al.l[2006h . During the first 3 s, the power of 
the injected outflow is Lj_o = PjTjfiVrj{hjTjfi — l)c^A5 = 
10^"'^ erg s~^, where AS* is the area of the injection sur- 
face, h{= 1 + tjc? + p/ pc?) is the specific enthalpy, and 
Vr is the radial component of the 3-vclocity. The density 
and pressure of the injected outflow are obtained by set- 
ting Fj^o, Ej, ij,o, 0^, andrmin- Wefixij^o = 10^^ erg s'^ 
which is high er than th a,t adopted in p revious studies 
(jZhang et all [20031 . [200I iMorsonv et al.ll2007n . The to- 
tal injected energy is several times 10^^ erg. Since the 



main purpose of this study is to sec the jet propaga- 
tion and expansion of the cocoon into the interstellar 
medium after the shock breakout, we adopt this power 
to obtain a rapid propagation of the jet in the progeni- 
tor. This fast propagation is necessary to be consistent 
with the fact that we neglect the self-gravity of the star. 
If the jet crosses the progenitor much faster than the 
typical hydrodynamic timescale in the system, the pro- 
genitor remains roughly unchanged during the complete 
jet propagation through it and, therefore, we do not need 
to care about the progenitor evolution during such short 
timescales. 

After the initial phase of constant kinetic power injec- 
tion, both the kinetic power and the injection Lorentz 
factor are linearly decreased according to the laws 
Lj(t) = max{Lj,o(4 - t),10''^crgs-i}, andFj(i) = 
max{Fj.o + 12— 4t, 1.01}, respectively, for 3s < t < 4s. In 
this period of decaying injection power, the specific en- 
ergy is kept fixed to the same value as it had at < = 0, and 
the density and the pressure are obtained from the other 
parameters (as in the constant injection power phase). 
After < = 4 s, the flow injection ceases. 

With the parametrization considered above, the rest- 
mass density (pj,o) during the constant power phase 
is 154 gcm""^. Since the rest- mass density of the pro- 
genitor around the inner computational boundary is 
lO^gcm"'^, the injected outflow is initially much lighter 
than medium in which it is injected. Thus it is expected 
that the jet propagation velocity across the progenitor 
is smaller than the speed of light, and one naturally ex- 
pects to generate relatively thick cocoons surrounding 
the beam of the jet. 

3. RESULTS 
3.1. Dynamics 

The dynamical evolution of our jet models can be split 
in two phases. The first one happens during the period in 
which the jet drills its way through the progenitor star. 
The second one shows up latter, when the jet breaks out 
of the stellar surface. The dynamics of our models during 
this two phases if roughl y similar to that outlined by 
some previ ous works (e.g .-IAlov et al.ll2000HZhang et all 
[2003, 2004: lMizuta et al. '2006: Morsonv et al.ll2007[ ) and, 
therefore, we limit ourselves here to provide a shallow 
description of the most salient features. 

Figure [H shows a snapshot of the evolution of the den- 
sity of the model IIE16N at t = 1.0 s, when the head 
of the jet is still in the progenitor. The left panel of 
the Fig. [H shows the Lorentz factor contour at the same 
time shown in Fig. [3] The jet is well collimatcd both in- 
side of the progenitor and as it travels through the ISM. 
The bow shock develops close to the head of the jet and 
rises the pressure and the temperature of the envelope 
region it sweeps up (in agreement with the findings of 
IMizuta et ani2004[ ). It takes about 3.2s for the jet to 
cross the progenitor, hence, the average propagation ve- 
locity is ~ 0.63c. The right panel of Fig. [4| shows the 
Lorentz factor of the model IIE16C at i = 0.7 s, show- 
ing that the head of the jet in model IIE16C propagates 
faster than model HE16N due to the lower density in 
model HE16C. 

Although the outflow has a finite initial opening angle, 
the the beam of the jet is almost parallel to the polar 
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HE16Nt=1.0 s 
Log10 (density) 




Fig. 3. — Snapshot (at t = 1.0 s) of the rest-mass density of model 
HE16N. The left and right panels are shown with different color 
scales in order to outline the most salient features of the cocoon 
and the cavity drilled by the jet (left), and to show more clearly 
the structure of the beam (right). Both, the vertical and the hori- 
zontal axis are scaled by Rq = Zq = 10^ cm. A strong bow shock 
surrounds the jet, the cocoon, and the shocked progenitor gas. A 
vigorous back flow from the head of the jet can be seen. Some vor- 
tices caused by Kelvin-Helmholtz instabilities are also indicated in 
the figure. 



axis. The kinetic energy of the beam is dissipated when 
it crosses the reverse shock (i.e, the Mach disc) at the 
head of the jet. After the beam plasma is decelerated at 
the Mach disc and its pressure is risen to a much higher 
value than in the beam, it expands and flows back in 
a thick cocoon. The high pressure of the cocoon is the 
responsible for the beam collimation during the initial 
phase of propagation inside of the progenitor star. Also 
during this early stage of the evolution, a strong backflow 
can be seen flanking the beam of the jet. Some vortices 
develop between the jet and the backflow caused by the 
growth of Kelvin-Helmholtz modes. An schemati c view 
of thi s process can be seen in Fig. 5a of Mizut a~et al.l 

dlooe'). 

The propagation of the jet inside the progenitor also 
drives a cavity limited by a shroud whose density and 
pressure is larger than in the cocoon. The shroud is 
swept up by a reverse shock that results from the inter- 
action between the cavity and the progenitor envelope. 
However, this reverse shock is not strong enough to rise 
the temperature above the threshold in which nuclear 
reactions can take place. 

After the jet breaks out the progenitor surface, it pro- 
ceeds to the ISM, which is assumed to be rarefied for 
model HE16N. In this phase, cocoon is almost freely re- 
leased into the ISM, because of the negligible pressure of 



HE16N HE16C 
t=1.0s t=0.7s 




Fig. 4. — Lorcntz factor contours of two models (left: model 
HE16N at t = 1.0s, and right HE16C at t= 0.7s). Both, the 
vertical and the horizontal axis are scaled by Ro = Zq = 10^ cm. 
The largest Lorentz factors are reached close to the polar axis where 
the first reconfinement shock appears. Some weak reconfinement 
shocks can also be seen further downstream. 

the external medium (Fig. In spite of the fact that 
the inertial confinement provided by the stellar progen- 
itor is lost in the ISM, by the time the jet reaches the 
surface of the star, the beam has accelerated to F > 40 
(see Sect. 13.1. ip and, thus, it has entered into a ballistic 
regime, where lateral expansion is strongly suppressed. 
Thereby, the jet remains well collimated as it propagates 
through the ISM, and the half-opening angle of the beam 
reaches only a few degrees. These collimation properties 
have been confirmed by means of numerical models with 
better resolution in the 0-direction (see Appendix). 

3.1.1. Acceleration to High Lorentz Factor 

In this section we focus on the dynamics of the accel- 
eration of the flow to large Lorentz factors F > 1 00, as 
reque sted by the the standard fireball model (e.g.. lPiraiil 
119991 ). Within such a theoretical model, an initial release 
of thermal energy is later converted into kinetic energy 
of the fiow, as it expands into a dilute and cold environ- 
ment. On the basis of this model, a number of numerical 
results have assessed that a jet, injected through a pre- 
established nozzle, is able to drill its way through a col- 
lapsar can reach an ultrarelativistic regime under likely 
inflow conditions (IZhang et al.ll2003l . 120041 : iMizuta et at] 
I2006t iMorsonv et all 120071 and this paper). We focus 
here on the details of the dynamical phase in which a 
kinematically mildly relativistic jet (Fj ~ 5) speeds up 
to the ultrarelativistic regime F > 100 converting its ini- 
tial (thermal) energy ej/c^ ^ 1 into kinetic energy. 
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HE16N t=4.0 s 
LoglO(density) Log10(r) 




R/Ro 



Fig. 5. — Density and Lorontz factor contours of the model 
HE16N at 4 = 4 s. Both, the vertical and the horizontal axis 
arc scaled by Rq = Zq = 10^ cm. By this time, the head of the 
jet has already broken out of the progenitor surface, and the bow 
shock begins to expand into the ISM. 



To drive the discussion, we take as a prototype case 
that of model HE16N. The acceleration process of other 
models is very similar. Figures [6^-[6l: show one dimen- 
sional profiles of density, pressure, and Lorentz factor 
along the polar axis, at different times. These profiles are 
quali tatively siniilar to the ones sh own in previous papers 
fe.g.. lAlov et al.ll2002l : iMorsonv et al. 2007,1. During the 
injection the Lorentz factor increases linearly (Fig. [G};), 
whereas both the density and pressure decrease as 
and (Fig. [S^, b), respectively. This is not unex- 
pected, since the fiuid expands radially (almost freely) 
in that region, where there are no shocks. Note that 
this result differs a bit from m odels where the genera - 
tion of the outflow is considered (|Alov et al.ll2000l . l2002f ) . 
If jet injection conditions are set through a nozzle at a 
certain distance to the center, the variability imprinted 
by the highly dynamical generation of the jet is erased. 
Clearly, this minimize s the number of internal shocks 
in the outflowing jet. lAlov et all ()200G[ ) show in their 
Fig. 2 that the outflow can accelerate to Lorentz fac- 
tors which are smaller than those attained in this work 
and in others where jet injection conditions are assumed. 
This is because of the modulation in the growth of the 
Lorentz factor imposed by the development of Kelvin- 
Helmholtz modes in the co urse of the very early outflow 
evolution (jAlov et al.ll2002f l. Nevertheless, models which 
consistently include the outflow generation (by thermal 
energy deposition), can still accelerate to ultrarelativistic 
Lorentz factors at large distance from the source. In such 
models the flow is kept hotter than in the present ones 




1e+08 1e+09 1e+10 1e+11 

radius r (cm) 

Fig. 6. — Profiles of the rest-mass density (a), of the pressure (b), 
of the bulk Lorentz factor (c), and of hT (d) along the polar axis 
at times t =1.0, 3.0, 3.5, 4.5 and 5.0s (see legends), corresponding 
to model HE16N. 



because of the occurrence of internal shocks already in- 
side of the progenitor. However, since the outfiow is opti- 
cally thick, the thermal energy is not r adiated away, but 
converted (later) into kinetic energ y JAlov et all [20001 : 
lAlovfc Mart3[200l : iGhisellini et al\\20m . 

One can notice that the maximum Lorentz factor is 
reached behind the first recoUimation shock in the jet 
(Figs, m and [Bt) coinciding with the location of the min- 
imum density in the beam of the jet (Figs. [3] and [5^). 
Inside of the progenitor star the jet reaches a maximum 
Lorentz factor F > 100 (see dashed line in Fig. [BJ: corre- 
sponding to t = 3 s) . We note that the quick acceleration 
to large Lorentz factor happening in the present models 
is a direct result of the fact that th e flow is injected int o 
a flnite opening angle. Differentlv. IMizuta et al.l lj2006f) . 
who injected the jet parallel to the polar axis of the pro- 
genitor, did not find such a fast acceleration. In the 
later case, the occurrence of recoUimation shocks much 
closer to the injection nozzle prevents the development 
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a freely expanding, unshocked jet extending so far away 
as in the models of this work. Indeed, such internal, rec- 
oUimation shocks are the responsible of the confinement 
of th e jet (Komissarov & Fallc '1997, 19981; iMizuta et aP 
12004 I2OOG: Mimic a ct, al. 2008 ). N onetheless, th e head 
propagation speed is very similar in lMizuta et al.l ()2006l ) 
as it is found here and, henc e, also the crossing time of 
the progenitor by the jet. In IMizuta et al.l (|2006( ) an ac- 
celeration of the jet to r > 100 at larger distances from 
the source happens (by thermal-to-kinetic energy con- 
version), but it still takes place inside of the progenitor 
star. This resembles more the acceleration process of out- 
flows generated self-consistently from the inner engine as 
pointed out above. 

The absolute maximum Lorcntz factor (Fmax ^ 130) 
is attained by the time we begin to reduce the luminos- 
ity of the inflow (at t = 3s, see § 12. 2p . By this time, 
the head of the jet has already broken up the surface of 
the progenitor. This numerical value is smaller than the 
maximum one which could be potentially reached, pro- 
vided the conditions of the injected outflow, according 
to which Fmax ~ Tj^o^j — 205, which is an indication of 
the maximum potentially reachable Lorentz factor, along 
a flow line, if all thermal energy were converted to ki- 
netic energy. After the energy deposition is switched 
off {t = 3 s), the region previously occupied by the un- 
shocked beam of the jet (up to the first rccoUimation 
shock of the beam; Fig. [6^) begins to shrink and to in- 
corporate mass from its lateral boundaries (c.f., rest-mass 
density profiles between times 3.5 s and 5.0 s in Fig. [6^ 
up to r ~ 8 X 10^ cm). This fact reduces progressively 
both the beam Lorentz factor (Fig. ^) and the ability 
of the beam fluid to reach large asymptotic values of F 
(Fig. [6jl). Thereby, if the jet is injected with a finite 
half-opening angle, the ultrarelativistic part of the out- 
flow (i.e., that where F > 100) may persist only if the 
activity of the engine does not cease. 

3.2. Angular Energy Distribution 

The evolution of our models has been followed until 
the head of the jet reaches the outermost radial compu- 
tational boundary at r,„ax ~ 10^^ cm. At such distances, 
the outflow is still optically thick and, therefore, radi- 
ation cannot escape freely. The jet and the surround- 
ing cocoon, past a transient phase, that happens after 
the jet breaks out of the stellar surface, enter into a 
quasi self-similar phase, where the properties of the out- 
flow roughly scale with the distance of propagation and 
become almost time independent. This fact allows us 
to extrapolate the properties of the jet from distances 
r ~ 10^"'^ cm to the typical distances where the after- 
glow will take place (namely, r ~ 10^'^ cm). Thus, wc 
can roughly estimate the angular energy distribution per 
unit solid angle {dE/d^l) that can be potentially emitted 
at the afterglow phase. 

In order to derive dE/dfl we have to make several 
assumptions. First, we assume that a fix fraction (the 
same everywhere) of the total energy (internal plus ki- 
netic) will be converted into electromagnetic radiation. 
Basically, this assumption is equivalent to state that 
the angular profile of the observed non-thermal radia- 
tion is simply a scaled version of the total energy an- 
gular profile. Certainly, this is a rough approximation, 
since the non-thermal radiation from 7-rays to radio fre- 



quencies will be produced by synchrotron (and, perhaps, 
inverse Compton) processes of particles accelerated at 
shocks (or, maybe, al ong the jet boundary layer; e.g., 
lAlov fc RezzoUal 120061 ) . Obviously, there are shocks of 
very different properties in the ultrarelativistic beam and 
in the cocoon and, thus, we may expect somewhat dif- 
ferent conversion efficiencies of the outflow energy into 
radiation in the beam and in the cocoon. Finally, we as- 
sume that the angular energy distribution is frozen-in by 
the time when the head of the jet reaches the outer com- 
putational boundary. As commented above, our models 
evolve almost self-similarly a bit after they break out of 
the stellar surface, and therefore, we expect only a minor 
time evolution of the angular profiles of dE/dCl. 

We point out that the pr ocedure we use to es timate 
dE/dD, differs from that of iMorsonv et all ()2007f ). who 
derived their dE/dVL profiles from the time integration 
of the energy flux trough a certain radius. Under the 
hypothesis of self-similar evolution, this is equivalent to 
integrate, along the radial direction, the energy density 
of our models (by the time they reach rmax) as follows. 



\ 1 — p--,. COS \o — Ul.) 

k=-Ng 1 = 1, ^ ^^1^ ^ 

^ ■" PikTik)ri'^^n, (1) 



X [pikhikTik 



1 - 

l-/3[,cos( 



where the subscripts i and k are associated to the spher- 
ical grid coordinates , and 9k , respectively, and and 
9 are the radial velocity in units of c and the observer's 
viewing angle (measured from the jet axis), respectively. 
The expression ([1]) includes the radiation contributions 
coming from reg ions outside of the line of sight (see 
iJanka et al.l 120061 ). The summation in the radial direc- 
tion runs from the surface of the progenitor, located at 
r = or, equivalently, z = z*, to the outermost bound- 
ary. The summation in the azimuthal angle runs from 
9 = -90° to 6* = 90° (note that due to the assumed 
axial symmetry, we can copy the computed data of the 
quadrant 0° < 6* < 90° to the quadrant -90° < 61 < 0°). 
In order to avoid accounting for subrelativistic regions, 
which will not contribute to the afterglow energetics, 
we exclude the contributions of numerical cells where 
Vr < 0.7c and ft.F < 4 in the expression 

The absolute value of the observed dE/dVl{9) along 
every radial direction forming an angle 9 with the po- 
lar axis depends, among other things, on two parameters 
whose exact value is not well constrained, neither by ob- 
servations nor by the present day theory. These are (i) 
the efficiency of energy conversion to radiation, and (ii) 
the total energy injected. Therefore, we will show only 
the angular profiles of dE/dil{9) normalized to the maxi- 
mum value dE / dn{9)\^_^^^ found for each model. Figure[7] 
shows the normalized angular energy distributions corre- 
sponding to models HE16C, HE16L and HE16N, which 
are prototypes of the types L, M and H, respectively. 
In the same figure we overplot fits to the normalized 
dE/dVl{9) profiles. The fitting function is a smoothly 

Note that numerical cells where hV < 4 have the potential to 
accelerate, at most, to F ~ 4. In actuality, the asymptotic Lorentz 
factor of such parcels of fluid, will be much smaller, since they will 
decelerate as they incorporate mass from the external medium. 
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0.01 - 



1e-04 




40.0 



Fig. 7. — Angular energy distribution of models, HE16C, HE16L, 
and HE16N, when the head of the jet reaches the outer computa- 
tional boundary. With lines we also show the fitting functions to 
SBPs (Eq. [2]l for each model. The solid, dashed, and dotted lines 
are the fitting functions of Eg. [2l corresponding to models HE16C, 
HE16L, and HE16N, respectively. The inset displays a zoom of 
the angular range 5° < 6 < 20°, where simple power law fitting 
functions are overplotted, along with the corresponding values of 
the power law indices. This angular region is dominated by the 
contribution of the expanding cocoon. 



broken power law (SBP) of the form 

_ ccin / n \ "h" 

F{e) ^ 2-i/"A 



l/r, 



(2) 



where A is the value of the function F at 9 = Oq, 9o the 
angular location of the break point between the prebreak 
and postbreak power-laws, whose slopes are ai and ah, 
respectively, and n is a numerical factor that controls the 
sharpness of the break. Note that the maximum value of 
F occurs at 



7max — Of) 



ah 

Oil 



l/[(a;-ah)n] 



(3) 



when aiUh < 0. Otherwise, if a; < and < 0, the 
function diverges as ^ 0. 

By inspection of Fig. [7l the angular energy distribu- 
tions are remarkably well fitted by the function of Eq. ([2]) 
in the interval < < 3.4°, i.e., in the angular region 
occupied by the beam of the jet. At smaller latitudes 
(5° 1^ S < 8°) the model data separates from the fit- 
ting function and presents systematically larger values 
than the latter. Indeed, the data in such an interval can 
be well fitted by a simple power law, with a slope in the 
range [—2.1, —2.5] (see the inset in Fig. [7]). The deviation 
from the SPB function in this angular range is due to the 
contribution of the expanding, mildly relativistic cocoon. 
This cocoon contribution shows up more clearly (the en- 
ergy distribution tends to flatten in the range 9 > 5°) if 
we lower the thresholds on the values of hT and w,. used to 
compute the angular energy distribution of each model. 
However, as we argued before, lowering these thresholds 
too much will catch up contributions from numerical cells 
whose asymptotic Lorentz factor is too small to account 
for typical afterglow. 

Low metallicity models 16TB, 16TC, and 160C, are 




0.01 - 



1e-04 



0.1 0.5 1.0 5.0 10.0 20.0 40.0 
angle (degree) 

Fig. 8.— Same as Fig.0 but for the models, 16TB, 16TC, and 
160C. 



also well fit by the function of Eq. ^ in roughly the 
same interval as the solar metallicity models (Fig. [8]). 
The values of the fit parameters are comparable to those 
of typc-H models, with which they share a very similar 
progenitor mass (~ ISM©). However, in these models, 
the cocoon contribution, which can be fit by a simple 
power-law with a slope ^ —2.6 between 5° < 6* < 8° 
(Fig. [8] inset) shows a faster decay of dE/dfl than type- 
H models for 9 > 8° . The reason for this difference is 
the much deeper density drop of low-metallicity models 
close to the star surface (Fig. [2|) compared with type-H 
models (Fig. [T]). The density of low-metallicity models 
in the region 3 x 10^°cm< r < 4 x 10^° cm is ^ 100 
times smaller than in type-H models. Hence, the beam 
of jets in such low-metallicity progenitors becomes much 
more ballistic than the corresponding beams of jets in 
the type-H group. Since more ballistic beams reduce the 
sideways expansion of their cocoons, this explains that 
the angular energy distribution in low-metallicity stars 
is more narrowly concentrated than in solar-mctallicity, 
typc-H progenitors. 

In order to show more clearly the existence of cor- 
relations between the properties of the progenitor star 
and the dE/dO. distribution, we show in Fig. [9] the de- 
pendence of the postbreak slope ah and on the stellar 
progenitor mass M . There exists a correlation between 
ah and M, such that the slope of lighter progenitors is 
steeper than that of heavier ones. There is roughly a lin- 
ear dependence of ah on M, which displays a relatively 
large dispersion. The reason for the dispersion being 
that for very similar values of the total progenitor mass, 
the rest-mass density radial profiles can be appreciably 
different (see, e.g., Figs. [T] and [2]). This is particularly 
true in heavy progenitor models (including type-H and 
low metallicity models). For the prebreak slope ai we 
find no obvious correlation with the progenitor mas, but 
in all the models considered here is very small (a/ ~ 0; 
Tab.[T|). 

We have also investigated the dependence of the slope 
ah and ai on the mass lo ss rate M assumed in models of 
iWooslev fc Heg«l ([20061 ) (see Tab. [J). Figure [TO] shows 
that there exist a good correlation between Uh and M, 
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Fig. 9. — Dependence of the index on the total mass of the 
progenitor. We identify the models by the last letter in the model 
name, e.g., the label 'A' stands for model HE16A. The labels 'TB', 
'TC, and 'OC stand for the models 16TB, 16TC, and 160C, re- 
spectively. 
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Fig. 11. — Equivalent isotropic angular energy distribution of 
the jets when the head of the jet reaches the outer computational 
boundary. The models shown are HE16C, HE16L, and HE16N, 
along with their corresponding fitting functions (with lines). 
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Fig. 10. — Dependence of the index on th e mass loss rate in 
units of the typical mass loss rate considered bv lWooslev fc Hegeil 
II2006I ). The naming convention is the same as in Fig.|9] 
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Fig. 12. — Same as Fig[9]but for the equivalent isotropic energy 
distribution. 



while ai seems to be independent of M. The ah — M cor- 
relation tells us that, models with a larger mass loss rate 
posses a steeper slopes. This is not surprising consider- 
ing the previously found correlation between and M , 
since the stellar progenitor mass is mostly determined by 
the amount of mass lost in the form of winds during the 
latest stages of its evolution. 

Wc have not found any other good correlation between 
the fit parameters (other than ah) and the gross prop- 
erties of the progenitors (radius, average density, total 
angular momentum, rotation period, mass of the iron 
core, etc.). 

We shall note that the angular distribution dE/dO, is 
not directly observable. Instead, the isotropic equivalent 
angular energy per solid angle dE / d^\iso can be detected. 
Figure [Tl] shows the equivalent isotropic angular energy 
distribution for prototype models belonging to the types 
L, M and H, normalized to the value of the distribu- 
tion at 6 = 0.25°. As can be seen from Fig. [TT] (see 
also Tab.[T|), the values of ai^iso are also negative for the 



dE / dD,\iso, while the corresponding values ai are close 
to zerofor dE/dQ (approximately, a/jso— cti — 1). This 
happens because of the small value of the solid angles 
close to the symmetry axis, which makes systematically 
larger the higher latitude values of dE/dn\iso than those 
ofdE/dQ. The values of ck/i. iso are systematically smaller 
than the respective ah values (roughly, ah.iso— ct/i — 1 in 
our standard resolution runs). Given the tight relations 
of the slopes in the dE/d^l and dE/d^l\iso distributions, 
it is not surprising to find that there exists a good corre- 
lation between ah,iso and M (Fig. [T^ which follows the 
same qualitative trend as the correlation between ah and 
M. 

To sum up our findings so far, we realize that the val- 
ues of the parameters of the fit function (Eql2|) are chiefiy 
correlated with the mass of the progenitor. In type- 
L progenitors, the velocity of propagation of the jet is 
larger than in more massive models, which results in jets 
developing more massive and hotter cocoons in type-M 
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and type-H models than in light progenitors^. This fact 
is evident by looking at Fig. 1131 where the pressure in 
the cocoon of the model HE16N is higher than that in 
the model HE16C in the course of the propagation of the 
jet up to the progenitor surface. The differences in the 
propagation velocity of the injected jets are set, to a large 
extend, by the differences in the average density between 
progenitors of different mass. The less massive progen- 
itors of our sample (type-L models) tend to have the 
smaller average densities (see, e.g.. Fig. [U where model 
HE16C displays a smaller density than model HE16N 
at every radial point). Consistently, jets propagating 
in type-L progenitors are faster as can be observed in 
Fig. [131 However, the time needed to reach the progen- 
itor surface presents a dependence with the progenitor 
mass with a much larger scatter (Fig. [TS)) . Particularly, 
models HE16A, HE16F, HE16G and HE16H, aU of which 
belong to the type-H group, display a progenitor crossing 
time comparable to that of the jets in the type-L group 
(although its average propagation velocity is comparable 
to that of models of type-H). That's the reason why these 
models appear as outliers in the correlations between the 
slope ah and the progenitor mass (Figs. [9] and with the 
mass- loss rate (Fig. I10|) . Since more massive progenitor 
stars yield hotter cocoons, once the cocoon is erupted 
through the stellar surface, it undergoes a larger lateral 
expansion (compare models HE16C and HE16N in the 
lower panels of Fig.[T]), i.e., the energy carried by the jet 
spreads towards lower latitudes. This explains why the 
energy per solid angle is more concentrated towards the 
axis when the jet crosses a low mass progenitor. 

If we assume that the detectability of an event, for 
an observer looking such event at a certain viewing 
angle (9), is pro portional to dE / dil\iso{d) (see, e.g., 
iJanka et al.ll200l . it turns out that 9 should be rather 
small (^2°; Fig. [T6|) to observe and event produced in 
collapsar progenitors (note that the flanks of dE/dn\iso 
distribution are quite steep and, thus, it is very unlikely 
to detect events which are not directly pointing towards 
the observer). Jets produced in type-L progenitors ex- 
hibit narrower observability profiles than those injected 
in more massive starts (c.f., compare the profiles of mod- 
els HE16C and HE16N in Fig.[ll|). Therefore, it is more 
unlikely to detect off-axis events produced in light pro- 
genitors than in more massive ones. Alternatively, we 
may state the the lower degree of collimation of relativis- 
tic jets in high-mass stars results in a higher probability 
of observing an event from a high-mass progenitor than 
a low-mass progenitor. We also find that the metallic- 
ity has little influence on the observability, because for 
similar progenitor masses, solar-metallicity models and 
low metallicity ones display almost identical observabil- 
ity angular profiles (Fig. fTB]). 

We also notice the very different observability profiles 
of jets produced in coUapsars and jets produced in rcm- 

* Since the injected mass flux is the same in aU models, jets 
which take longer to reach the stellar surface (i.e., jets with small 
propagation speed) are more massive and, hence their cocoons are 
also heavier. On the other hand, the slower jets develop stronger 
reverse shocks than the faster ones, because the injected fluxes of 
mass, momentum and energy are the same in all our models. Thus, 
the jet matter is compressed more in the former than in the latter 
case, which explains why the matter in the cocoon is hotter in jets 
with smaller propagation speeds. 
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Fig. 13. — Pressure contours of models HE16C (left: at f = 2.0 s) 
and HE16N (right: at t = 3.2 s), when the head of the jet reaches 
progenitor surface in each model, i.e., at shock break out. Both, the 
vertical and the horizontal axis are scaled by Rq = Zq = 10^ cm. 
The radii of the progenitors are 4.73 X 10^" cm and 6.13X 10^'^ cm for 
the models HE16C and HE16N, respectively. At the same distance 
Z/Z{i from the injection nozzle, the pressure in both the jet and 
in the cocoon is higher for model HE16N than for model HE16C. 
This fact explains the larger lateral expansion after shock break 
out in model HE16N than in model HE16C. 
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Fig. 14. — Dependence on the progenitor mass of the average 
propagation velocity of the jet in the star. The average propagation 
velocity is deflned as the ratio between the progenitor radius r» and 
the time it takes the head of the jet to propagate up to r, . The 
labeling convention is the same as in Fig. [§] 



Relativistic Jet Eruption from CoUapsars 



11 







*N 












•TC 






• 1 ■ 


B-K 
J *C 




tA 


6 8 

— Dependence 


10 12 
progenitor mass (M^^^) 

on the progenitor 


14 16 
mass of the time 



Fig. 15. 

ployed by the jet to reach the progenitor surface, 
convention is the same as in Fig. |9] 
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Fig. 16.— Dctectability of the models HE16C, HE16L, HE16N 
and 16TB. 



nants of mergers of compact objects. In the latter case, a 
non-negligible observability is obtained in the viewing an- 
gle interval 2°<6'<13°, because of the larger half-opening 
angles of the outflows originated in progenitors of short 
GRBs (see Fig. 2 of Janka ct al. 2006). Another relevant 
difference is that collapsar-jcts show a dctectability with 
a fast rise close to 6* = 0, followed by a shallower decay 
beyond the most probable detection angle ^ 0.5°. In 
contrast, most jets produced in merger remnants tend to 
decay very abruptly beyond the most likely observing an- 
gle, and have an observability rise much more moderate 
than that of coUapsar-jets. The reason for this difference 
is the much shallower decay of the isotropic equivalent 
energy at the flanks of the beam in coUapsar-jets than in 
jets from merger remnants. 

4. CONCLUSIONS AND DISCUSSION 

In this paper we explore the relation of the dynamical 
properties of ultrarelativistic jets generated in coUapsars 
with the properties of the progenitor stars in which such 
jets propagate. We have particularly focused on the cor- 
relations that exist between the angular proflle of the 
energy per solid angle (as seen in the afterglow phase) 



and the properties of the progenitors. Along the way, we 
have pointed out which is the relevance of the fact that 
our numerical models are set up with a finite injection 
half-opening angle. 

Using a non-zero injection angle affects the way in 
which the conversion of internal-to-kinetic energy takes 
place. If the flow is injected parallel to the polar axis, 
the development of reconfinement shocks happens closer 
to the injection nozzle than if the flow is injected radi- 
ally within a cone of finite half-opening angle. When 
the rccoUimation shock occurs far away from the nozzle, 
the unshocked beam fiow accelerates along a larger dis- 
tance in a rarefaction that precedes such shock. Thus, 
the beam reaches there larger Lorentz factors than if 
the jet is injected parallel to the polar axis. The dis- 
sipation in cross shocks acts by simply recycling part of 
the kinetic energy of the outfiow into thermal energy. 
The thermal energy is not lost, since the jet propagation 
is roughly adiabatic inside of the progenitor (radiation 
losses are negligible there). Instead, this thermal energy 
can be further converted into kinetic energy at larger 
distances. This process may happen several times be- 
fore the outflow becomes transparent and radiation can 
freely escape. This explains why, in spite of the differ- 
ences in the beam dynamics, all mo dels (independent o f 
the injection half-opening angle; c.f.. lMizuta et al]|2006f) 
develop a roughly similar propagation speed and, by the 
time they reach the head of the jet, the gross properties 
of the outflow are similar. 

We have estimated the angular distribution of energy 
per solid angle in the afterglow phase by extrapolation 
of the state of our models when they reach a distance 
of ~ 10^^ cm. This extrapolation relies on the fact that 
the jets develop a rough self-similar behavior soon af- 
ter they emerge from the progenitor surface. Our re- 
sults show that the equivalent isotropic energy per solid 
angle dE/dVl\iso is only pa rtly consistent with that of 
iLazzati fc BegelmanI ()2005[ ). However, the results in 
this paper do appear to be consiste nt with previous nu - 



mer ical simulations such as those of I Zhang et al 
and IMorsonv et all (|2007f ). ILazzati fc Begelmaii 



( 2004 ) 



20051 ) 



obtain that the angular energy distribution displays a 
relatively flat core which is flanked by a region where 
dE/<KV\iso oc . Our results show that the core of the 
distribution (close to 6* = 0°) is not flat (but decays as 
6^^) and that the energy per solid angle decays much 
faster than 9^^ (it does it as 0'', with a value of k < —3.6 
depending on the mass of the progenitor; see below). 
We can fit the dE / dri\iso-da,t& with SEP functions up to 
9 < 3.4°. At smaller latitudes, a simple power-law with 
a slope close to —2.6 fits better the data. In this region 
(4.5° < 9 < 8°), the cocoon contribution is the dominant 
(at smaller values of 9, the beam of the jet dominates the 
energetics), and we find there the best c onsistency with 
the model of ILazzati fc BegelmanI (j2005D . 

We have correlated the properties of the angular 
dE/dQ\iso distribution of the jet with the fundamental 
parameters of the progenitor star in which the jet has 
propagated. We find that the shape of the distribution 
is mostly influenced by the mass of the progenitor. When 
the mass of the progenitor is small (A/ 5Mq), because 
of the occurrence of large mass losses due to winds in 
the latest stages of the star's evolution, then dE/dn\iso 
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decays faster with than if the mass of the progenitor is 
large (M ^ 15Mq). We find that the reason for this be- 
havior is that the average density of the progenitors tends 
to grow (approximately) with the mass. This means that 
the average jet propagation speed inside of the star is 
smaller, the larger is the mass. A smaller jet propaga- 
tion speed results into thicker and hotter cocoons, since 
we fix the same mass, momentum and energy fluxes at 
the nozzle for all models. Also the beam of the jet in the 
more massive progenitors is wider. This is the reason 
why low mass progenitors develop narrower dE/dCl\iso 
profiles than high mass ones. The difference in the col- 
limation of the energy distributions resulting from low 
and high mass progenitors has a direct influence on the 
number of observed events. We expect to see more events 
produced in heavy progenitors than in low mass ones. 

One could question whether the correlation that we 
have found between the mass of the progenitor and the 
width of the dE/dn\[so is an artifact of our numerical set 
up. We are fixing the luminosity of the jet to be the same 
independent of the progenitor mass. However, progeni- 
tors of different mass may develop central engines which 
release different power. A good proxy of the power of the 
central engine is the mass of the iron core of the progen- 
itor. One may expect that the collapse of more massive 
iron cores results into larger central compact objects. If 
the power released by the central engine is dominated 
by the size of the central compact object, then models 
with more massive iron cores could release a larger power 
than models with low mass cores. Then, how do we jus- 
tify our numerical assumption that the power injected in 
the jet is roughly independent of the progenitor's mass?. 
In support of o ur point we argue that, first, according 
to iWooslev fc H egcr ( 20061) data, there is not a one to 
one correlation between the mass of the iron core and 
the mass of the progenitor and, second, the iron mass 
varies by less than a 30% in all the models considered 
here, while the total mass can be different by a factor of 
3. Thus, within the simplifications we do in our models, 
the assumption of a common luminosity independent of 
the mass of the progenitor is justified. 

Irrespective of these two arguments, if heavier progen- 
itors would result into more luminous central engines, 
we shall point out that this trend will also result into 
wider j ets and cocoons. This is a result early pointed 
out bv lAlov et al.l (|200G[ ). were it was shown that, tak- 
ing the same progenitor, but increasing the luminosity of 
the central engine by a factor of 10, results into thicker 
cocoons than in cases in which the luminosity is more 
moderate. The reason being that the release of a large 
power triggers large amplitude Kelvin-Helmholtz insta- 
bilities at the basis of the outflow, which transfer a large 
fraction of the momentum of the beam to a thicker shear 



layer between the beam and the cocoon. Effectively, this 
process widens the cross sectional area of the beam, re- 
ducing its propagation speed and, consistently inflating 
larger cocoons. 

We have also found, that comparing progenitors of sim- 
ilar mass, the metallicity of the star has a small impact 
on the extrapolated dE/dfl\iso profiles. The coUimation 
of the jet is similar regardless of the stellar metallicity. 
However, the cocoon is more narrowly collimated in low- 
metallicity stars, because of the large density drop close 
to their surfaces. This reduced density makes the jets 
in low- metallicity stars much more ballistic, once they 
break out the stellar surface, than in solar-metallicity 
progenitors. Unfortunately, this difference might not be 
observable, since it happens in regions where the en- 
ergy per solid angle is much smaller than at the jet 
core (unless orphan afterglows could be detected; see, 
e.g., iTotani fc Panaitescul 120021 : iRossi. Perna fc Daignd 
[20081) . 

Finally, we have found significant differences between 
the dE/dil\iso profiles of coUapsar-jets and those of jets 
produced in merger remnants (i.e., between the angu- 
lar energy distribution of jets associated to long and to 
short GRBs). Collapsar-jets are more narrowly colli- 
mated than jets of merger remnants, and the decay of 
dE/dD,\iso beyond the central flat core is much steeper 
in the latter than in the former. This intrinsic difference 
manifest itself as a larger chance of detectability of jets 
from merger remnants at viewing angles up to ~ 12°, 
while, on the other hand, collapsar-jets could hardly be 
seen off-axis. 
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APPENDIX 

ON THE CHOICE OF NUMERICAL RESOLUTION 

The main purpose of this Appendix is to justify the choice of numerical resolution that we have used in the main 
body of this paper. We have picked up three models, HE16C, HE16L and HE16N (representative of the models of the 
respective types-L, -M and H) and performed a resolution study by progressively increasing the numerical resolution. 
Our standard models have a working resolution N,. x Ng = 1000 x 180 zones. The standard computational grid is 
uniform in both logr- and ^-coordinates. In addition to the standard resolution, two higher resolutions have been 
considered. 

First, we compute models HE16C-M, HE16L-M, and HE16N-M (Tab. [T|), which have an intermediate resolution of 
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HE16C HE16C-M HE16C-H 

t=2.0s t=1.9s t=2.1 s 




Fig. 17. — Contours of the energy density for the same model run at three different resolutions, corresponding to models HE16C, HE16C- 
M, and HE16C-H. The jet in all cases is about to break the surface of the star (_R» = 4.75 X 10 cm). Since the jet propagates at sHghtly 
different speed depending on the resolution, the three snapshots correspond to slightly different evolutionary times. Both, the vertical and 
the horizontal axis are scaled by Ro = Zq = 10^ cm. 



Nr X Ng ^ 1500 X 180 zones. In these models the radial grid is uniform in logr and the smallest radial grid spacing 
is Ar = 7.5 x 10^ cm, i.e., 3/4 times smaller than that of our standard resolution cases. The polar grid possesses a 
uniform region close to the symmetry axis (0° < 6* < 40°), where A0 = 1/3°, followed by a uniformly spaced region in 
log6' (30° < 6 < 90°). The reason to consider two different regions in the ^-spacing is that all our jet models develop 
cocoons whose angular extension is 6* < 30°, i.e., all the dynamics develops in a wedge covered by a finer mesh in our 
computational grid. Thus, an increased resolution in the abovementioned wedge surrounding the polar axis yields an 
effective increase of the numerical resolution in the whole computational grid. 

Even higher resolution models (HE16C-H, HE16L-H, and HE16N-H) have also been computed. In this case, the grid 
consists of Nr = 2000 zones uniformly spaced in logr, with a minimum radial grid spacing Ar = 5 x 10^ cm, i.e., one 
half of the width of smallest radial spacing of the standard resolution case. We take in this case Ng = 180 zones, also 
split in two regions: a uniform grid in the interval 0° < 6* < 30° {A9 = 0.25°), followed by a uniformly spaced region 
in log0 (30° <e < 90°). 

Going to even larger resolutions in the ^-direction increases the total computational time up to prohibitive limits 
due to the increased number of time steps associated to the fulfillment of the Courant condition. 

We note that the jet dynamics is rather independent of the resolution, and also the gross morphological features are 
converged at the standard resolution, though, of course, finer details show up both in the cocoon and in the beam (see 
Fig. [17]). Therefore, when we take radial averages to compute the angular energy distribution of our models (Eq. [1]), 
the differences are relatively small (see below) , and our standard resolution models can be considered to be sufficiently 
resolved to account for such global energetic properties. The time-scales to cross the progenitor and/or to reach the 
outer computational boundary (at r = 10^^ cm) are slightly different from each other. 

Figure [18] shows the isotropic equivalent angular energy distributions of models computed with three different 
resolutions. In all cases, a SBP function (Eq. [2]) fits properly the data up to 6* < 3.4°. At larger latitude, the 
contribution of the expanding cocoon component dominates. Since each distribution is normalized independently to 
its absolute maximum, the distributions corresponding to different resolutions do not overlap, but they show the same 
shape. The fitting parameters for the higher resolution cases arc also listed in Tab. [TJ There, we can see that the same 
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angle (degree) 

Fig. 18. — Comparison of the angular energy distributions of prototype models (HE16C (top), HE16L (middle), and HE16L (bottom)) 
run with different computational resolutions. The solid, dashed, and dotted lines show the fitting functions of Eq. |2] corresponding to 
models computed with the standard, the middle, and the highest resolution, respectively. 



correlations found in models with the standard resolution are reproduced at higher resolutions, namely, the correlation 
between ah or a/j^sowith the progenitor mass (Sect. 13. 2p . However, at higher resolution, we could guess additional 
correlations which are not obvious in models run with the standard resolution. Particularly, the parameters 9q and 
^o,iso become smaller as the progenitor mass increases. This trend is however, an artifact of the models chosen as 
prototypes of e ach mass type. T hey are such, that the radius of the progenitor roughly grows with the mass. As 
pointed out bv lAlov et alTpOOOf) . progenitors with larger radii provide a larger inertial confinement which prevents 
the lateral expansion of the jet. Hence, for the models chosen in this resolution study, a larger progenitor mass yields 
better collimatcd jets. Taking the whole sample of models, but run at higher resolution, the correlations between 
or ^o.iso would not exist. 
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TABLE 1 

Properties of the pre-supernova models taken from iWoosley HegerJ i|2006I'i for our study. Columns 1-4, list the model 

NAME, THE MASS LOSS RATE (iN THE SAME UNITS AS IWOOS LEY HeGER 2OO60, THE TOTAL MASS, AND THE RADIUS AT PRE-SUPERNOVA 

STAGE. Columns 5-9, display the best parameters of the SBP (Eq.Q used to fit the angular energy distribution per solid 
ANGLE (Eq.[TJ. Finally, columns 10-14, show the best fit parameters for the angular distribution of equivalent isotropic 

ENERGY PER SOLID ANGLE. THE LAST SIX ROWS CORRESPOND TO SOME PROTOTYPE MODELS WHICH HAVE BEEN RUN WITH HIGHER 
NUMERICAL RESOLUTION. MODELS WHOSE NAME ENDS WITH "-M" AND "-H" HAVE BEEN RUN WITH NUMERICAL RESOLUTIONS OF 1500 X 180 

AND 2000 X 180 GRID POINTS, RESPECTIVELY (SEE ApPENDIx). 
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